DEV ROC betaseries

load packages

library(tidyverse)
library(lme4)
library(lmerTest)
library(knitr)
library(cowplot)
library(caret)
library(ROCR)

define aesthetics

algorithm = c("#006989", "#FEC601", "#F43C13", "#00A5CF", "#00A878")
instruction = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
craving = wesanderson::wes_palette("Darjeeling1", 3, "continuous")
rating = c("#00A08A", "#F2AD00", "#F98400", "#FF0000")
dc_bw = plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank())

load tidied data

source("load_data.R")

check ratings

Check if wrong buttons were used (i.e., not 5-8)

• DEV001 = code normally
• DEV011 = code normally
• DEV016 = code normally
• DEV017 = exclude; can’t tell if they’re missed ratings or incorrect placement of fingers
• DEV019 = exclude; can’t tell if they’re missed ratings or incorrect placement of fingers
• DEV020 = code normally
• DEV022 = code normally
• DEV028 = code normally
• DEV032 = incorrect placement of fingers; recode runs 1-2 (LOOK INTO WTP)
• DEV037 = exclude; technical error?
• DEV054 = exclude; technical error?
• DEV060 = code normally; task ended early
• DEV061 = code normally; task ended early
• DEV063 = code normally; task ended early
• DEV069 = incorrect placement of fingers in run1
• DEV075 = code normally
• DEV082 = code normally
• DEV083 = code normally

NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES

subs = data.all %>%
  group_by(subjectID, wave, run, rating) %>%
  summarize(n = n()) %>%
  spread(rating, n) %>%
  mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
  filter(messed == "yes") %>% 
  ungroup() %>% 
  select(subjectID, wave) %>% 
  unique()

data.all %>%
  group_by(subjectID, run, rating) %>%
  summarize(n = n()) %>%
  spread(rating, n) %>%
  mutate(messed = ifelse(is.na(`5`) & !is.na(`<NA>`), "yes", NA)) %>%
  filter(subjectID %in% subs$subjectID)

recode and exclude

Recoding
* DEV069: recode runs1

NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES

data.ex = data.all %>%
  mutate(rating = ifelse(subjectID == "DEV069" & run == "run1", rating - 1, rating),
         rating = ifelse(subjectID == "DEV069" & run == "run1" & is.na(rating), 8, rating),
         rating = rating - 4) %>%
  group_by(subjectID, wave) %>%
  arrange(subjectID, run) %>%
  mutate(trial = row_number())

load mean intensity values

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave1/"
file_pattern = "DEV[0-9]{3}_meanIntensity.txt"
file_list = list.files(file_dir, pattern = file_pattern)

intensities = data.frame()

for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "meanIntensity" = V3) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    mutate(beta = as.integer(beta),
                           wave = 1), error = function(e) message(file))
  intensities = rbind(intensities, temp)
  rm(temp)
}

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave2/"
file_list = list.files(file_dir, pattern = file_pattern)
for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "meanIntensity" = V3) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    mutate(beta = as.integer(beta),
                           wave = 2), error = function(e) message(file))
  intensities = rbind(intensities, temp)
  rm(temp)
}

load dot products

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave1/"
file_pattern = "DEV[0-9]{3}_dotProducts.txt"
file_list = list.files(file_dir, pattern = file_pattern)

dots = data.frame()

for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "map" = V3,
                           "dotProduct" = V4) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    extract(map, "algorithm", "(.*)_.*.nii") %>%
                    mutate(beta = as.integer(beta),
                            wave = 1), error = function(e) message(file))
  dots = rbind(dots, temp)
  rm(temp)
}

file_dir = "~/Documents/code/sanlab/DEV_scripts/fMRI/betaseries/ROC/dotProducts_ROC_wave2/"
file_list = list.files(file_dir, pattern = file_pattern)

for (file in file_list) {
  temp = tryCatch(read.table(file.path(file_dir,file), fill = TRUE) %>%
                    rename("subjectID" = V1,
                           "map" = V3,
                           "dotProduct" = V4) %>%
                    extract(V2, "beta", "beta_([0-9]{4}).nii") %>%
                    extract(map, "algorithm", "(.*)_.*.nii") %>%
                    mutate(beta = as.integer(beta),
                            wave = 2), error = function(e) message(file))
  dots = rbind(dots, temp)
  rm(temp)
}

join intensities and dots

• recode trials with extreme intensities as NA

dots.merged = dots %>%
  left_join(., intensities, by = c("subjectID", "wave", "beta")) %>%
  group_by(subjectID, wave, algorithm) %>%
  mutate(rownum = row_number())

# plot original
dots.merged %>%
  filter(algorithm == "craving_regulation") %>%
  ggplot(aes(1, meanIntensity)) +
    geom_boxplot()

# assess extreme values and exclude when calculating SDs
dots.merged %>%
  filter(algorithm == "craving_regulation") %>%
  arrange(meanIntensity)

dots.merged %>%
  filter(algorithm == "craving_regulation") %>%
  arrange(-meanIntensity)

# recode outliers as NA
dots.merged = dots.merged %>%
  ungroup() %>%
  mutate(meanIntensity = ifelse(meanIntensity > 1 | meanIntensity < -1, NA, meanIntensity),
         median = median(meanIntensity, na.rm = TRUE),
         sd3 = 3*sd(meanIntensity, na.rm = TRUE),
         outlier = ifelse(meanIntensity > median + sd3 | meanIntensity < median - sd3, "yes", "no"),
         dotProduct = ifelse(outlier == "yes", NA, dotProduct))
  
# plot after
dots.merged %>%
  filter(algorithm == "craving_regulation") %>%
  ggplot(aes(1, meanIntensity)) +
    geom_boxplot()

recode subs

• DEV022 = run4 has 8 trials
• DEV037 = ???
• DEV048 = run4 missing
• DEV060 = run1 has 19 trials; couldn’t estimate run1 trial 19, run3 trial 20
• DEV061 = run3 has 19 trials; couldn’t estimate run3 trial 19
• DEV063 = run2 has 11 trials
• DEV081 = run2 missing (run1 was run twice)
• DEV082 = run2 has 15 trials; couldn’t estimate run1 trial 19, run1 trial 20

NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES

trial.numbers = data.frame(subjectID = c(rep("DEV060", 79), rep("DEV061", 79), rep("DEV063", 71), rep("DEV081", 80), rep("DEV082", 75)),
                           wave = 1,
                           rownum = c(1:79, 1:79, 1:71, 1:80, 1:75),
                           trial = c(1:19, 21:80, 1:59, 61:80, 1:31, 41:80, 1:20, 41:80, 21:40, 1:35, 41:80))

dots.check = dots.merged %>%
  group_by(subjectID, wave, algorithm) %>%
  mutate(rownum = row_number()) %>%
  left_join(., trial.numbers, by = c("subjectID", "wave", "rownum")) %>%
  mutate(trial = ifelse(is.na(trial), rownum, trial),
         dotProduct = ifelse(subjectID == "DEV060" & wave == 1 & trial %in% 19:20, NA,
                      ifelse(subjectID == "DEV061" & wave == 1 & trial == 59, NA,
                      ifelse(subjectID == "DEV082" & wave == 1 & trial %in% 19:20, NA, dotProduct)))) %>%
  select(-rownum) #%>%
  #left_join(., striping, by = c("subjectID", "beta")) %>%
  #mutate(dotProduct = ifelse(!is.na(striping), NA, dotProduct))

exclude outliers and standardize

• standardize within subject and algorithm

dots.ex = dots.check %>%
  group_by(subjectID, algorithm) %>% # standardize within sub and algorithm
  mutate(dotSTD = scale(dotProduct, center = FALSE)) 

merge data and exclude subs

Exclusions

• MRI motion and data quality exclusions: DEV001, DEV020, DEV032, DEV047, DEV055, DEV064, DEV066
• Button box exclusions: DEV017, DEV019, DEV037, DEV054
• Run exclusions: DEV029 (run3), DEV037 (run1), DEV042 (run4), DEV067 (run4)

Other
* select only craved trials

NB! NEED TO GO THROUGH NEW PARTICPANTS/WAVES

data_all = left_join(dots.ex, data.ex, by = c("subjectID", "wave", "trial")) %>%
  filter(!(subjectID %in% c("DEV001","DEV020","DEV032","DEV047","DEV055","DEV064","DEV066", "DEV017", "DEV019", "DEV037", "DEV054") & wave == 1)) %>%
  filter(!(subjectID == "DEV029" & wave == 1 & run == "run3") & !(subjectID == "DEV037" & wave == 1 & run == "run1") & !(subjectID == "DEV042" & wave == 1 & run == "run4") & !(subjectID == "DEV067" & wave == 1 & run == "run4")) %>%
  ungroup() %>%
  mutate(algorithm = gsub("_signature", "", algorithm),
         wave = as.character(wave))

data = data_all %>%
  filter(craving == "craved")

summarize

n participants by wave

data %>%
  select(subjectID, wave) %>%
  unique() %>%
  group_by(wave) %>%
  summarize(n = n())

n trials

data %>%
  filter(algorithm == "craving_regulation") %>%
  group_by(subjectID) %>%
  summarize(n = n()) %>%
  arrange(n)

roc

craving regulation signature

# roc curve
perf.df = data %>%
  filter(algorithm == "craving_regulation") %>%
  filter(!is.na(dotSTD) & !is.na(craving)) %>%
  mutate(instruction = ifelse(instruction == "regulate", 1, 0)) %>%
  group_by(wave) %>%
  do({
    wave = .$wave
    pred = prediction(.$dotSTD, .$instruction, label.ordering = NULL)
    perf = performance(pred, measure = "tpr", x.measure = "fpr")
    data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
  })

ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
  geom_line() +
  scale_color_manual(values = algorithm) +
  xlab("false positive rate") +
  ylab("true positive rate")

roc = data %>%
  filter(algorithm == "craving_regulation") %>%
  select(subjectID, trial, instruction, rating, dotSTD) %>%
  mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
         guess.rating = ifelse(dotSTD > 0, "low", "high"),
         rating.bin = ifelse(rating >= 3, "high", "low"),
         instruction = as.factor(instruction),
         guess.instruction = as.factor(guess.instruction),
         rating.bin = as.factor(rating.bin),
         guess.rating = as.factor(guess.rating))

confusionMatrix(roc$guess.instruction, roc$instruction)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction look regulate
##   look     5551     2292
##   regulate 3490     6707
##                                                
##                Accuracy : 0.6795               
##                  95% CI : (0.6726, 0.6863)     
##     No Information Rate : 0.5012               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.3592               
##                                                
##  Mcnemar's Test P-Value : < 0.00000000000000022
##                                                
##             Sensitivity : 0.6140               
##             Specificity : 0.7453               
##          Pos Pred Value : 0.7078               
##          Neg Pred Value : 0.6577               
##              Prevalence : 0.5012               
##          Detection Rate : 0.3077               
##    Detection Prevalence : 0.4348               
##       Balanced Accuracy : 0.6796               
##                                                
##        'Positive' Class : look                 
## 

#confusionMatrix(roc$guess.rating, roc$rating.bin)

craving signature

crave regulate v craved look

# roc curve
perf.df = data %>%
  filter(algorithm == "craving") %>%
  filter(!is.na(dotSTD) & !is.na(craving)) %>%
  mutate(instruction = ifelse(instruction == "regulate", 1, 0)) %>%
  group_by(wave) %>%
  do({
    wave = .$wave
    pred = prediction(.$dotSTD, .$instruction, label.ordering = NULL)
    perf = performance(pred, measure = "tpr", x.measure = "fpr")
    data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
  })

ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
  geom_line() +
  scale_color_manual(values = algorithm) +
  xlab("false positive rate") +
  ylab("true positive rate")

roc = data %>%
  filter(algorithm == "craving") %>%
  select(subjectID, trial, instruction, rating, dotSTD) %>%
  mutate(guess.instruction = ifelse(dotSTD > 0, "regulate", "look"),
         guess.rating = ifelse(dotSTD > 0, "low", "high"),
         rating.bin = ifelse(rating >= 3, "high", "low"),
         instruction = as.factor(instruction),
         guess.instruction = as.factor(guess.instruction),
         rating.bin = as.factor(rating.bin),
         guess.rating = as.factor(guess.rating))

confusionMatrix(roc$guess.instruction, roc$instruction)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction look regulate
##   look     3629     4178
##   regulate 5412     4821
##                                              
##                Accuracy : 0.4684             
##                  95% CI : (0.4611, 0.4757)   
##     No Information Rate : 0.5012             
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : -0.0629            
##                                              
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.4014             
##             Specificity : 0.5357             
##          Pos Pred Value : 0.4648             
##          Neg Pred Value : 0.4711             
##              Prevalence : 0.5012             
##          Detection Rate : 0.2012             
##    Detection Prevalence : 0.4328             
##       Balanced Accuracy : 0.4686             
##                                              
##        'Positive' Class : look               
## 

#confusionMatrix(roc$guess.rating, roc$rating.bin)

craving v neutral

# roc curve
perf.df = data_all %>%
  filter(algorithm == "craving") %>%
  filter(!is.na(dotSTD)) %>%
  filter(craving %in% c("neutral", "craved")) %>%
  mutate(craving = ifelse(craving == "craved", 1, 0)) %>%
  group_by(wave) %>%
  do({
    wave = .$wave
    pred = prediction(.$dotSTD, .$craving, label.ordering = NULL)
    perf = performance(pred, measure = "tpr", x.measure = "fpr")
    data.frame(cut=perf@alpha.values[[1]],fpr=perf@x.values[[1]],tpr=perf@y.values[[1]])
  })

ggplot(perf.df, aes(fpr, tpr, color = as.factor(wave))) +
  geom_line() +
  scale_color_manual(values = algorithm) +
  xlab("false positive rate") +
  ylab("true positive rate")

roc = data_all %>%
  filter(algorithm == "craving") %>%
  filter(craving %in% c("neutral", "craved")) %>%
  select(subjectID, trial, craving, rating, dotSTD) %>%
  mutate(guess.craving = ifelse(dotSTD > 0, "craved", "neutral"),
         guess.rating = ifelse(dotSTD > 0, "low", "high"),
         rating.bin = ifelse(rating >= 3, "high", "low"),
         craving = as.factor(craving),
         guess.craving = as.factor(guess.craving),
         rating.bin = as.factor(rating.bin),
         guess.rating = as.factor(guess.rating))

confusionMatrix(roc$guess.craving, roc$craving)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction craved neutral
##    craved   10233    3897
##    neutral   7807    5143
##                                              
##                Accuracy : 0.5678             
##                  95% CI : (0.5619, 0.5737)   
##     No Information Rate : 0.6662             
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.1229             
##                                              
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.5672             
##             Specificity : 0.5689             
##          Pos Pred Value : 0.7242             
##          Neg Pred Value : 0.3971             
##              Prevalence : 0.6662             
##          Detection Rate : 0.3779             
##    Detection Prevalence : 0.5218             
##       Balanced Accuracy : 0.5681             
##                                              
##        'Positive' Class : craved             
## 

#confusionMatrix(roc$guess.rating, roc$rating.bin)

visualize raw data

instruction

across waves

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(algorithm, dotSTD, fill = instruction)) +
    stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
    scale_fill_manual(name = "", values = instruction) + 
    labs(y = "standardized pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(instruction, dotSTD)) +
    stat_summary(aes(group = subjectID), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = 1), fun.y = mean, geom = "line", size = .75) +
    stat_summary(aes(color = instruction), fun.data = "mean_cl_boot",  geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = instruction) + 
    labs(y = "standardized pattern expression value\n", x = "") + 
    dc_bw

by waves

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(instruction, dotSTD, fill = wave)) +
    stat_summary(fun.y = mean, geom = "bar", position = position_dodge(width = 0.95)) +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.95), width = 0) +
    facet_grid(~algorithm) +
    scale_fill_manual(name = "wave", values = algorithm) + 
    labs(y = "standardized pattern expression value\n", x = "") + 
    dc_bw

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(wave, dotSTD, color = instruction)) +
    stat_summary(aes(group = interaction(subjectID, instruction)), fun.y = mean, geom = "line", alpha = .1, size = .5) +
    stat_summary(aes(group = instruction), fun.y = mean, geom = "line", size = .75) +
    stat_summary(fun.data = "mean_cl_boot",  geom = "pointrange", width = 0, size = .75) + 
    facet_grid(~algorithm) +
    scale_color_manual(name = "", values = instruction) + 
    labs(y = "standardized pattern expression value\n", x = "") + 
    dc_bw

rating

across waves

data %>%
  ggplot(aes(dotSTD, rating, color = algorithm, fill = algorithm)) +
  geom_jitter(alpha = .1, height = .1) +
  geom_smooth(method = "lm") +
  scale_color_manual(name = "signature", values = algorithm) + 
  scale_fill_manual(name = "signature", values = algorithm) + 
  labs(y = "craving rating\n", x = "\nstandardized pattern expression value") + 
  dc_bw

by wave

data %>%
  ggplot(aes(dotSTD, rating, color = wave, fill = wave)) +
  geom_jitter(alpha = .1, height = .1) +
  geom_smooth(method = "lm") +
  facet_grid(~algorithm) +
  scale_color_manual(name = "wave", values = algorithm) + 
  scale_fill_manual(name = "wave", values = algorithm) + 
  labs(y = "craving rating\n", x = "\nstandardized pattern expression value") + 
  dc_bw

rating and instruction

across waves

data %>%
  ggplot(aes(dotSTD, rating, color = instruction, fill = instruction)) +
  geom_jitter(alpha = .1, height = .1) +
  geom_smooth(method = "lm") +
  facet_grid(~algorithm) +
  scale_color_manual(name = "instruction", values = instruction) + 
  scale_fill_manual(name = "instruction", values = instruction) + 
  labs(y = "craving rating\n", x = "\nstandardized pattern expression value") + 
  dc_bw

by wave

data %>%
  ggplot(aes(dotSTD, rating, color = wave, fill = wave)) +
  geom_jitter(alpha = .05, height = .1) +
  geom_smooth(method = "lm") +
  facet_grid(instruction~algorithm) +
  scale_color_manual(name = "wave", values = algorithm) + 
  scale_fill_manual(name = "wave", values = algorithm) + 
  labs(y = "craving rating\n", x = "\nstandardized pattern expression value") + 
  dc_bw

RT

across waves

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(dotSTD, rt, color = instruction, fill = instruction)) +
  geom_jitter(alpha = .05, height = .1) +
  geom_smooth(method = "lm", alpha = .2) + 
  facet_grid(~algorithm) + 
  scale_color_manual(name = "", values = instruction) + 
  scale_fill_manual(name = "", values = instruction) + 
  labs(x = "\nstandardized pattern expression value", y = "reaction time (seconds)\n") + 
  dc_bw

data %>%
  filter(!is.na(rating)) %>%
  filter(algorithm == "craving_regulation") %>%
  ggplot(aes(rating, rt, color = instruction, fill = instruction)) +
  geom_jitter(alpha = .05, height = .1) +
  geom_smooth(method = "lm", alpha = .2) + 
  facet_grid(~algorithm) + 
  scale_color_manual(name = "", values = instruction) + 
  scale_fill_manual(name = "", values = instruction) + 
  labs(x = "\ncraving rating", y = "reaction time (seconds)\n") + 
  dc_bw

by wave

data %>%
  filter(!is.na(rating)) %>%
  ggplot(aes(dotSTD, rt, color = wave, fill = wave)) +
  geom_jitter(alpha = .05, height = .1) +
  geom_smooth(method = "lm", alpha = .2) + 
  facet_grid(instruction~algorithm) + 
  scale_color_manual(name = "", values = algorithm) + 
  scale_fill_manual(name = "", values = algorithm) + 
  labs(x = "\nstandardized pattern expression value", y = "reaction time (seconds)\n") + 
  dc_bw

data %>%
  filter(!is.na(rating)) %>%
    filter(algorithm == "craving_regulation") %>%

  ggplot(aes(rating, rt, color = wave, fill = wave)) +
  geom_jitter(alpha = .05, height = .1) +
  geom_smooth(method = "lm", alpha = .2) + 
  facet_grid(~instruction) + 
  scale_color_manual(name = "", values = algorithm) + 
  scale_fill_manual(name = "", values = algorithm) + 
  labs(x = "\ncraving rating", y = "reaction time (seconds)\n") + 
  dc_bw

run models

craving regulation signature

pattern expression ~ wave * instruction

When looking, people showed increases in expression of the craving regulation signature (i.e., their brains looks more like they’re regulating) from wave 1 to 2.

When regulating, people showed decreases in expression of the craving regulation signature (i.e., their brains look less like they’re regulating) from wave 1 to 2. This may mean that it’s gotten easier and/or they’re exprience less craving and therefore need to regulate less.

mod_pev_wave = lmer(dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID),
                    data = filter(data, algorithm == "craving_regulation"))

ggeffects::ggpredict(mod_pev_wave, terms = c("wave", "instruction")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(group = group)) +
  scale_color_manual(name = "instruction", values = instruction) + 
  scale_fill_manual(name = "instruction", values = instruction) + 
  labs(x = "\nwave", y = "predicted pattern expression value\n") + 
  dc_bw

summary(mod_pev_wave)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID)
##    Data: filter(data, algorithm == "craving_regulation")
## 
## REML criterion at convergence: 46494.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1474 -0.6296  0.0024  0.6300  4.4455 
## 
## Random effects:
##  Groups    Name                      Variance Std.Dev. Corr             
##  subjectID (Intercept)               0.07721  0.2779                    
##            wave2                     0.05818  0.2412   -0.37            
##            instructionregulate       0.10015  0.3165   -0.04  0.03      
##            wave2:instructionregulate 0.11528  0.3395    0.12 -0.38 -0.38
##  Residual                            0.72303  0.8503                    
## Number of obs: 18040, groups:  subjectID, 255
## 
## Fixed effects:
##                            Estimate Std. Error        df t value
## (Intercept)                -0.27366    0.02138 252.53513 -12.801
## wave2                       0.08587    0.02445 229.88478   3.513
## instructionregulate         0.89880    0.02645 247.41990  33.983
## wave2:instructionregulate  -0.15568    0.03440 221.59634  -4.525
##                                       Pr(>|t|)    
## (Intercept)               < 0.0000000000000002 ***
## wave2                                 0.000534 ***
## instructionregulate       < 0.0000000000000002 ***
## wave2:instructionregulate           0.00000985 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wave2  instrc
## wave2       -0.480              
## instrctnrgl -0.290  0.245       
## wv2:nstrctn  0.264 -0.554 -0.512

craving rating ~ wave * pattern expression * instruction

Cravings decreased from wave 1 to 2; the effect was stronger when people were looking compared to when they were regulating.

The more the brain looks like a person is regulating, the lower their cravings; this relationships is somewhat stronger (i.e., more negative) when people are looking compared to regulating.

mod_craving = lmer(rating ~ wave * instruction * dotSTD + (1 + wave * instruction | subjectID),
                    data = filter(data, algorithm == "craving_regulation"))

ggeffects::ggpredict(mod_craving, terms = c("dotSTD", "wave", "instruction")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
  geom_line(aes(group = group)) +
  facet_grid(~facet) +
  scale_color_manual(name = "wave", values = algorithm) + 
  scale_fill_manual(name = "wave", values = algorithm) + 
  labs(x = "\npredicted pattern expression value", y = "predicted craving rating\n") + 
  dc_bw

ggeffects::ggpredict(mod_craving, terms = c("wave", "dotSTD[-1,0,1]", "instruction")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
  geom_line(aes(group = group)) +
  facet_grid(~facet) +
  scale_color_manual(name = "pattern expression", values = algorithm) + 
  scale_fill_manual(name = "pattern expression", values = algorithm) + 
  labs(x = "\nwave", y = "predicted craving rating\n") + 
  dc_bw

summary(mod_craving)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rating ~ wave * instruction * dotSTD + (1 + wave * instruction |  
##     subjectID)
##    Data: filter(data, algorithm == "craving_regulation")
## 
## REML criterion at convergence: 37063.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8913 -0.6678  0.0148  0.6644  4.0521 
## 
## Random effects:
##  Groups    Name                      Variance Std.Dev. Corr             
##  subjectID (Intercept)               0.2142   0.4628                    
##            wave2                     0.3191   0.5649   -0.17            
##            instructionregulate       0.2775   0.5268   -0.45 -0.06      
##            wave2:instructionregulate 0.1752   0.4186   -0.11 -0.58 -0.28
##  Residual                            0.4505   0.6712                    
## Number of obs: 17148, groups:  subjectID, 253
## 
## Fixed effects:
##                                     Estimate  Std. Error          df t value
## (Intercept)                          3.25803     0.03132   253.04648 104.027
## wave2                               -0.60967     0.04175   216.05262 -14.604
## instructionregulate                 -0.92645     0.03723   272.05521 -24.885
## dotSTD                              -0.03429     0.01213 16012.60486  -2.827
## wave2:instructionregulate            0.10576     0.03732   268.45316   2.834
## wave2:dotSTD                        -0.01295     0.01764 16177.48242  -0.734
## instructionregulate:dotSTD           0.02846     0.01677 15853.36189   1.697
## wave2:instructionregulate:dotSTD     0.02790     0.02411 14963.52335   1.157
##                                              Pr(>|t|)    
## (Intercept)                      < 0.0000000000000002 ***
## wave2                            < 0.0000000000000002 ***
## instructionregulate              < 0.0000000000000002 ***
## dotSTD                                        0.00471 ** 
## wave2:instructionregulate                     0.00494 ** 
## wave2:dotSTD                                  0.46287    
## instructionregulate:dotSTD                    0.08970 .  
## wave2:instructionregulate:dotSTD              0.24715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wave2  instrc dotSTD wv2:ns w2:STD in:STD
## wave2       -0.235                                          
## instrctnrgl -0.475  0.036                                   
## dotSTD       0.108 -0.082 -0.093                            
## wv2:nstrctn  0.034 -0.580 -0.382  0.094                     
## wave2:dtSTD -0.074  0.099  0.064 -0.683 -0.115              
## instrct:STD -0.077  0.058 -0.070 -0.718  0.069  0.489       
## wv2:nst:STD  0.053 -0.071  0.048  0.493 -0.103 -0.723 -0.685

craving signature

pattern expression ~ wave * instruction

Pattern expression of the craving signature decreased from wave 1 to 2 (i.e., peoples’ brains looked less like they were craving).

Regulation is associated with lower craving signature expression across waves.

mod_pev_wave = lmer(dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID),
                    data = filter(data, algorithm == "craving"))

ggeffects::ggpredict(mod_pev_wave, terms = c("wave", "instruction")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(group = group)) +
  scale_color_manual(name = "instruction", values = instruction) + 
  scale_fill_manual(name = "instruction", values = instruction) + 
  labs(x = "\nwave", y = "predicted pattern expression value\n") + 
  dc_bw

summary(mod_pev_wave)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: dotSTD ~ wave * instruction + (1 + wave * instruction | subjectID)
##    Data: filter(data, algorithm == "craving")
## 
## REML criterion at convergence: 49649.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.2765 -0.6227  0.0061  0.6223  6.7929 
## 
## Random effects:
##  Groups    Name                      Variance Std.Dev. Corr             
##  subjectID (Intercept)               0.08926  0.2988                    
##            wave2                     0.08692  0.2948   -0.44            
##            instructionregulate       0.04248  0.2061    0.15 -0.05      
##            wave2:instructionregulate 0.04861  0.2205    0.22 -0.57 -0.15
##  Residual                            0.87165  0.9336                    
## Number of obs: 18040, groups:  subjectID, 255
## 
## Fixed effects:
##                            Estimate Std. Error        df t value
## (Intercept)                 0.23606    0.02315 249.22640  10.199
## wave2                      -0.05843    0.02808 235.86125  -2.081
## instructionregulate        -0.13192    0.02299 249.39700  -5.737
## wave2:instructionregulate  -0.01078    0.03185 236.89089  -0.338
##                                       Pr(>|t|)    
## (Intercept)               < 0.0000000000000002 ***
## wave2                                   0.0386 *  
## instructionregulate               0.0000000278 ***
## wave2:instructionregulate               0.7353    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wave2  instrc
## wave2       -0.520              
## instrctnrgl -0.271  0.259       
## wv2:nstrctn  0.322 -0.620 -0.531

craving rating ~ wave * pattern expression * instruction

The more the brain looks like a person is craving, the higher their cravings (though this relationship was weak). This relationships is somewhat stronger when people are regulating compared to looking.

The relationship between craving signature pattern expression and craving was somewhat stronger at wave 2 compared to wave 1.

mod_craving = lmer(rating ~ wave * instruction * dotSTD + (1 + wave * instruction | subjectID),
                    data = filter(data, algorithm == "craving"))

ggeffects::ggpredict(mod_craving, terms = c("dotSTD", "wave", "instruction")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted, color = group, fill = group)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5, color = NA) +
  geom_line(aes(group = group)) +
  facet_grid(~facet) +
  scale_color_manual(name = "wave", values = algorithm) + 
  scale_fill_manual(name = "wave", values = algorithm) + 
  labs(x = "\npredicted pattern expression value", y = "predicted craving rating\n") + 
  dc_bw

summary(mod_craving)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rating ~ wave * instruction * dotSTD + (1 + wave * instruction |  
##     subjectID)
##    Data: filter(data, algorithm == "craving")
## 
## REML criterion at convergence: 37023.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9327 -0.6625  0.0156  0.6661  3.9978 
## 
## Random effects:
##  Groups    Name                      Variance Std.Dev. Corr             
##  subjectID (Intercept)               0.2147   0.4634                    
##            wave2                     0.3244   0.5696   -0.17            
##            instructionregulate       0.2793   0.5285   -0.45 -0.05      
##            wave2:instructionregulate 0.1797   0.4239   -0.09 -0.59 -0.29
##  Residual                            0.4492   0.6702                    
## Number of obs: 17148, groups:  subjectID, 253
## 
## Fixed effects:
##                                     Estimate  Std. Error          df t value
## (Intercept)                          3.26349     0.03128   250.60545 104.342
## wave2                               -0.61426     0.04196   215.24030 -14.639
## instructionregulate                 -0.94102     0.03655   250.48483 -25.743
## dotSTD                               0.01721     0.01127 16055.19923   1.527
## wave2:instructionregulate            0.12236     0.03612   229.32164   3.388
## wave2:dotSTD                         0.02940     0.01586 16328.59675   1.854
## instructionregulate:dotSTD           0.03003     0.01578 15837.07175   1.903
## wave2:instructionregulate:dotSTD    -0.02092     0.02197 15395.51824  -0.952
##                                              Pr(>|t|)    
## (Intercept)                      < 0.0000000000000002 ***
## wave2                            < 0.0000000000000002 ***
## instructionregulate              < 0.0000000000000002 ***
## dotSTD                                       0.126901    
## wave2:instructionregulate                    0.000829 ***
## wave2:dotSTD                                 0.063818 .  
## instructionregulate:dotSTD                   0.057023 .  
## wave2:instructionregulate:dotSTD             0.341120    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) wave2  instrc dotSTD wv2:ns w2:STD in:STD
## wave2       -0.235                                          
## instrctnrgl -0.483  0.040                                   
## dotSTD      -0.083  0.064  0.070                            
## wv2:nstrctn  0.040 -0.608 -0.368 -0.074                     
## wave2:dtSTD  0.058 -0.080 -0.049 -0.708  0.092              
## instrct:STD  0.058 -0.045 -0.074 -0.708  0.076  0.499       
## wv2:nst:STD -0.041  0.056  0.051  0.503 -0.087 -0.713 -0.709

within-person mediation

prep data

# source functions
library(boot)
source("indirectMLM.R")

# create craving signature dataframe
data_med_craving = data %>%
  filter(algorithm == "craving") %>%
  mutate(wave = ifelse(wave == "1", 0,
                ifelse(wave == "2", 1, NA)),
         dotSTD = dotSTD[,]) %>%
  select(subjectID, wave, trial, instruction, rating, dotSTD) %>%
  data.frame() %>%
  na.omit()

# create craving regulation signature dataframe
data_med_regulation = data %>%
  filter(algorithm == "craving_regulation") %>%
  mutate(wave = ifelse(wave == "1", 0,
                ifelse(wave == "2", 1, NA)),
         dotSTD = dotSTD[,]) %>%
  select(subjectID, wave, trial, instruction, rating, dotSTD) %>%
  data.frame() %>%
  na.omit()

# define variables
x_var = "wave"
y_var = "rating"
m_var = "dotSTD"

craving signature

all trials

model_name = "mediation_craving"
data_med = data_med_craving

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.001, 0.004]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.005 [-0.007, 0]
## Biased Estimate of Within-subjects Indirect Effect: -0.005 [-0.008, -0.002]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.545 [-0.571, -0.517]
## Biased Total Effect of X on Y (c path): -0.545 [-0.571, -0.516]
## Bias in Total Effect: 0 [0, 0.003]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.541 [-0.567, -0.513]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.064 [-0.09, -0.033]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.082 [0.067, 0.095]

look trials

model_name = "mediation_craving_look"
data_med = filter(data_med_craving, instruction == "look")

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.003]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.002 [-0.006, 0.001]
## Biased Estimate of Within-subjects Indirect Effect: -0.002 [-0.004, 0]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.609 [-0.637, -0.576]
## Biased Total Effect of X on Y (c path): -0.611 [-0.639, -0.579]
## Bias in Total Effect: 0.002 [0, 0.005]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.607 [-0.636, -0.574]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.053 [-0.092, -0.014]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.032 [0.014, 0.05]

regulate trials

model_name = "mediation_craving_regulate"
data_med = filter(data_med_craving, instruction == "regulate")

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.001]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.004 [-0.009, -0.002]
## Biased Estimate of Within-subjects Indirect Effect: -0.004 [-0.006, -0.002]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.494 [-0.522, -0.468]
## Biased Total Effect of X on Y (c path): -0.495 [-0.523, -0.469]
## Bias in Total Effect: 0.001 [0, 0.004]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.49 [-0.518, -0.462]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.073 [-0.118, -0.034]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.053 [0.038, 0.067]

craving signature

all trials

model_name = "mediation_regulation"
data_med = data_med_regulation

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.005, 0.002]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.002 [-0.008, 0.005]
## Biased Estimate of Within-subjects Indirect Effect: -0.002 [-0.005, 0.005]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.005]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.551 [-0.568, -0.531]
## Biased Total Effect of X on Y (c path): -0.545 [-0.561, -0.524]
## Bias in Total Effect: 0.006 [0.002, 0.007]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.549 [-0.564, -0.533]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.008 [-0.02, 0.022]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): -0.234 [-0.243, -0.227]

look trials

model_name = "mediation_regulation_look"
data_med = filter(data_med_regulation, instruction == "look")

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: -0.001 [-0.005, 0.002]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: -0.004 [-0.009, -0.001]
## Biased Estimate of Within-subjects Indirect Effect: -0.003 [-0.006, -0.002]
## Bias in Within-subjects Indirect Effect: 0.001 [0, 0.005]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.61 [-0.637, -0.572]
## Biased Total Effect of X on Y (c path): -0.611 [-0.637, -0.573]
## Bias in Total Effect: 0.001 [0, 0.004]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.606 [-0.634, -0.567]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.086 [0.053, 0.124]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): -0.041 [-0.061, -0.022]

regulate trials

model_name = "mediation_regulation_regulate"
data_med = filter(data_med_regulation, instruction == "regulate")

if (file.exists(sprintf("models/model_%s.RDS", model_name))) {
  assign(get("model_name"), readRDS(sprintf("models/model_%s.RDS", model_name)))
} else {
  assign(get("model_name"), boot(data = data_med, statistic = indirect.mlm, R = 500,
                                 y = y_var, x = x_var, mediator = m_var, group.id = "subjectID",
                                 between.m = F, uncentered.x = F))
  saveRDS(eval(parse(text = model_name)), sprintf("models/model_%s.RDS", model_name))
}

indirect.mlm.summary(get(model_name))

## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0 [-0.004, 0.003]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: 0 [-0.004, 0.003]
## Biased Estimate of Within-subjects Indirect Effect: 0 [-0.001, 0.001]
## Bias in Within-subjects Indirect Effect: 0 [0, 0.004]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: -0.495 [-0.521, -0.466]
## Biased Total Effect of X on Y (c path): -0.495 [-0.522, -0.467]
## Bias in Total Effect: 0 [0, 0.004]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): -0.494 [-0.52, -0.467]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): -0.062 [-0.097, -0.023]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.003 [-0.012, 0.021]